
% read out parameters to access them with their name
load ./markov2/output/markov2_results.mat

NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;

%%%%%%%%%% Don't change the above code %%%%%%%%%%%%
global theta dbar  L gL alpha rbest rho sigma ns n betas d alphas fspace_smolyak ilinear ilog_rL delta beta



% private eqm
r_private = rbest;
n  = ncountries;
ns = nsectors;
ilog_rL = 1; % use exp(rL)/(1+exp(rL))

dij = dbar;
d=dij*ones(n); %dij is (1+iceberg) cost producing in j selling to i
d=d-(dij-1)*eye(n);


for i=1:n
    eval(['L(',num2str(i),',1)=L',num2str(i),';'])
    for j=1:ns
        eval(['alpha(',num2str(i),',',num2str(j),')=alpha',num2str(i),num2str(j),';'])
    end
end

for j=1:ns
    eval(['betas(1,',num2str(j),')=betas',num2str(j),';'])
end

%% solve for the steady state of gov assuming dG2/dT=0!!!
% This is not the true steady state, it's only for us to get initial guess
% --- steady state ---


options = optimoptions('fsolve', 'TolFun', 1e-12);
options.Display = 'off';

xguess(1:2) = 1.1;
xguess(3) = 0.2;

[xval,fval_sol,eflag_ss] = fsolve(@(xval) solve_1s_2c_gov_ss_nogL(xval),xguess,options); %iter


[Tss,xss,Pss,pi_ss, fval_ss] = eval_1s_2c_gov_ss_nogL(xval);

xval_ss= xval;
wss = xval(1:n);
extax_ss = xval(3);
taf_ss = 0;
rss = rbest;

log(Pss(2))-(-1/theta*log(Tss(1)*(wss(1)*(1+extax_ss)*dbar)^(-theta)+Tss(2)*(wss(2))^(-theta)))


% construct M
for i=1:n
    GTss(i)= log(xss(i)^(-sigma)*Pss(i)^(sigma-1)*wss(i)/alpha(i));
    Mss(i) = wss(i)/alpha(i);
end

uc = xss(1)^(-sigma);

for i=2:n
    for j=1:ns
        for ii=1:n
            for jj=1:ns
                eval(['g0=gcoef0_',num2str(i),num2str(j),num2str(ii),num2str(jj),';']);
                dGT(i,ii) = g0*exp(GTss(i))/Tss(ii);
            end
        end
    end
end
               
gammaT_ss = wss(2)/alpha(2)*(1+theta)/theta*pi_ss(2,1)*uc*extax_ss/(1+extax_ss);
mu = gammaT_ss*alpha(2);

%% write to our variable names in dynare
for j=1:ns
    eval(['wedge',num2str(j),'=0;'])
end
rL = rbest*L;
Lp_ss=(1-rbest)*L;
ilog=0;

if (ilog==1)
fprintf(1,'w1 = %20.10e;\n',log(wss(1)));
fprintf(1,'w2 = %20.10e;\n',log(wss(2)));
fprintf(1,'T1 = %20.10e;\n',log(Tss(1)));
fprintf(1,'T2 = %20.10e;\n',log(Tss(2)));
fprintf(1,'pi12 = %20.10e;\n',pi_ss(1,2));
fprintf(1,'pi21 = %20.10e;\n',pi_ss(2,1));
fprintf(1,'x1 = %20.10e;\n',log(xss(1)));
fprintf(1,'x2 = %20.10e;\n',log(xss(2)));
fprintf(1,'r1 = %20.10e;\n',log(rbest));
fprintf(1,'r2 = %20.10e;\n',log(rbest));
fprintf(1,'taux = %20.10e;\n',extax_ss);
fprintf(1,'pi11 = %20.10e;\n',1-pi_ss(1,2));
fprintf(1,'pi22 = %20.10e;\n',1-pi_ss(2,1));
fprintf(1,'P2 = %20.10e;\n',log(Pss(2)));
fprintf(1,'M2 = %20.10e;\n',Mss(2));
fprintf(1,'M1 = %20.10e;\n',Mss(1));
else
    fprintf(1,'w1 = %20.10e;\n',wss(1));
    fprintf(1,'w2 = %20.10e;\n',wss(2));
    fprintf(1,'T1 = %20.10e;\n',Tss(1));
    fprintf(1,'T2 = %20.10e;\n',Tss(2));
    fprintf(1,'pi12 = %20.10e;\n',pi_ss(1,2));
    fprintf(1,'pi21 = %20.10e;\n',pi_ss(2,1));
    fprintf(1,'x1 = %20.10e;\n',xss(1));
    fprintf(1,'x2 = %20.10e;\n',xss(2));
    fprintf(1,'r1 = %20.10e;\n',rbest);
    fprintf(1,'r2 = %20.10e;\n',rbest);
    fprintf(1,'taux = %20.10e;\n',extax_ss);
    fprintf(1,'pi11 = %20.10e;\n',1-pi_ss(1,2));
    fprintf(1,'pi22 = %20.10e;\n',1-pi_ss(2,1));
    fprintf(1,'P2 = %20.10e;\n',Pss(2));
    fprintf(1,'M2 = %20.10e;\n',Mss(2));
    fprintf(1,'M1 = %20.10e;\n',Mss(1));
end

